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We performed dynamic simulations of spheres with short-range attractive interactions for many values of interaction strength 
and range. Fast crystallization occurs in a localized region of this parameter space, but the character of crystallization pathways 
is not uniform within this region. Pathways range from one-step, in which a crystal nucleates directly from a gas, to two-step, 
in which substantial liquid-like clusters form and only subsequently become crystalline. Crystallization can fail because of slow 
nucleation from either gas or liquid, or because of dynamic arrest caused by strong interactions. Arrested states are characterized 
by the formation of networks of face-sharing tetrahedra that can be detected by a local common neighbor analysis. 


Colloidal crystallization is of considerable interest because 
of the value of colloidal assemblies to technol^ypS and the 
value of colloidal dispersions as model systems!^. Colloidal 
particles can be made of controlled size, interaction strength, 
and interaction range, making them useful models for explor¬ 
ing the thermodynamic and kinetic factors that lead to the as¬ 
sembly of equilibrium and nonequilibrium condensed states of 
matter^. 

A defining feature of many colloidal suspensions is that their 
interactions can be short ranged compared to the nm to /im size 
of the colloidal particles. For example, van der Waals inter- 
acti on, depl etion interactions® and DNA base-pairing interac- 
tionsEQStni ^sed to promote colloidal crystallization typically 
act over a range of distances small compared to the particle 
size. As a result, colloidal suspensions can exhibit phase be¬ 
havior and assembly kinetics not typically seen in atomic or 
molecular systemst^HHI. 

A minimal and well-studied model of a colloidal disper¬ 
sion is a collection of spher es of radius Rq interacting via the 
‘square-weir potential^^^^^^"^ 

( oo rij<2RQ 

U{rij)=l -£ 2/?o<r,y <2(l+A)/?o (1) 

[ 0 r,-,->2(l+A)/?o, 

where rij is the distance between the centers of particles i and 
j. In this model the solvent in which particles are dispersed 
is not represented explicitly. Square-well spheres display the 
well -known phase behavior of mrticles with isotropic attrac- 
tionsEUmi^ summarized in Fig.^ When the interaction range 
A is large (Fig.[^(a)), the square-well system exhibits sequen¬ 
tial phase transitions from gas to liquid to crystal as temperature 
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(i.e. the combination k^T/e) decreases. (Here “gas” refers to 
a dilute suspension of colloids and “liquid” to a concentrated 
suspension.) As the interaction range decreases, the gas-liquid 
coexistence curve decreases in temperature (Fig. [^(b)), even¬ 
tually becoming metastable with respect to fiuid-crystal coex¬ 
istence (Fig.[^(c)). 

Many studies have shown that such metastable liquid-gas 
phase separation can play a cruc ial ro le in colloidal crystal- 
lizat ion. Free-en ergy calculationslSini and dynamic simula- 
short-range attractive spheres show that 
the metastable liquid can promote a two-step crystallization 
pathway in which colloids coalesce into liquid droplets from 
which crystals nucleate. Two-step pathways have been ob¬ 
served experimentally in colloidal particles confined in two di- 
mensioml 23 [ 28 |^ DNA-tethered nanoparticles®, and in pro- 
teinsEHEI]^ as well as in simulations of DNA-tethered nanopar¬ 
ticles® 

Colloidal liquid-gas phase separation also plays an impor¬ 
tant role in the formation of (physical) gels at deep supercool¬ 
ing. Gels are nonequilibrium, disordered networks of parti¬ 
cles with solid-like mechanical properties that result from their 
percolating structures. Gelation occurs because strong inter¬ 
particle bonding causes particles within aggregates to rearrange 
too slowly to allow equilibration on observed timescales. In 
the deeply supercooled spinodal regime, rapid liquid-gas phase 
separation can cause the formation of extended, non-compact 
colloid al agg regates which fail to relax into compact colloidal 
droplets®^. Many experimental studies of polymeric colloidal 
particles possessing depletion attractio ns ha ve found gelation 
to occur in preference to crystallizafiun®Sl (potentially exac¬ 
erbated by effects of polydispersity®®). Microscopic analy¬ 
sis of the colloid-colloid interaction networks formed during 
gelation shows gelation to be characterized by certain (over¬ 
lapping) locally-favored motifs: gels in long-range repulsive 
colloids consist of networks of face-sharing tetrahedra (max- 
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Figure 1 Phase diagrams of hard spheres of radius Rq with attractive square-well interactions of range 2(1 + in the plane of density 
(hard-core packing fraction 0) and temperature (scaled interaction strength k^T/e). Interaction ranges are (a) A = 0.33, (b) A = 0.25, and (c) 
A = 0.18. Blue circles represent densities of coexisting fluids, calculated using Gibbs ensemble Monte Carlo simulations, and blue curves are 
fits to the coexistence curve for systems in the Ising universality class. Black circles represent densities of fluids coexisting with 
face-centered-cubic (fee) crystals, calculated using direct coexistence simulations. Crystal densities, to the right of these plots, are not shown. 


imally bonded clusters of four particles)EISl, while gels of 
short-range-attractive spheres consist of networks of maximally 
bonded clusters of various sizes® 

Recently, several authors have investigated the generality 
of crystallization and gelation mechanisms by characterizing 
colloid dynamics across broad sections of parameter space. 
Macfarlane et al. showed in experiments that DNA-linked 
nanoparticle crystallization occurs for each nanoparticle size 
only within a limited range of DNA lengths. Short lengths re¬ 
sulted in effective interaction ranges smaller than the nanoparti¬ 
cles’ polydispersity, disfavoring the costal thermodynamically, 
while large lengths inhibited kinetics®. 

Several authors have performed Monte Carlo, molecular dy¬ 
namics, or Brownian dynamics simulations of spheres with 
short-range attractive interactions to investigate how assembly 
mechanism and product depend on interaction strength and/or 
concentration. Multiple studies have shown that crystallization 
occurs via two-step nucleation for a window of temperatures 
below the metastable liquid-gas tran sition, with deeper tem¬ 
perature quenches leading to gelation®2Illl^ Extending these 
studies to multiple concentrations, Fortini et al. showed that 
crystallization coincides with the metastable liquid-gas transi¬ 
tion temperature at concentrations below the liquid-gas criti¬ 
cal concentration but occurs also at higher temperatures at su¬ 
percritical concentrations®. Performing single-particle Monte 
Carlo simulations of square-well spheres at a single packing 
fraction of (j) = 0.04, Klotsa and Jack found an exception to 
the two-step rule: at a temperature near the metastable liquid- 
gas transition, they found that crystallization can proceed via a 
one-step pathway without significant formation of amorphous 


clusters®. 

Complementing and extending these studies, we describe in 
this paper the self-assembly dynamics of square-well-attractive 
spheres over a broad spectrum of interaction strengths and 
ranges. We simulated sphere dynamics using the virtual-move 
Monte Carlo algorithml^SESl parameterized so that colloidal 
clusters diffuse at rates agreeing with Stokes’ law. Consis¬ 
tent with the previously mentioned studies at fixed interaction 
range, we find that efficient crystallization occurs in a local¬ 
ized region of parameter space, with a high-temperature bound¬ 
ary associated with the metastable liquid-gas transition. How¬ 
ever, we find that the character of crystallization pathway varies 
within the region of efficient crystallization. Near the high- 
temperature boundary, crystallization proceeds along a one- 
step pathway via nucleation from the gas. Further below this 
boundary, crystallization proceeds along a two-step pathway 
via the formation of liquid-like clusters from which crystals 
subsequently nucleate. We find that poor crystallization at low 
temperature is characterized by the formation of networks of 
face-sharing tetrahedra that can be detected by a local common 
neighbor analysis®. 

It is important to note that the square-well model we have 
studied neglects features of real colloidal particles that may 
lead to complexity beyond that discussed here. The pairwise 
nature of the square-well interaction cannot capture collective 
properties of counterions, depletants, or polymer coats that me¬ 
diate multi-body interactions between colloidal particles. For 
example, counterion entropy can favor gelation over crystal¬ 
lization in a way that cannot be modeled at a pairwise level®, 
and the entropy of mobile linkers in the dilute-linker limit can 
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Figure 2 (a) Illustration of the common neighbor analysis for a 
bonded pair of spheres with a 423 common neighbor environment: 
the bonded pair shares four common neighbors, there are two bonds 
among those common neighbors, and three common neighbors 
participate in those two bonds, (b-g) Low-energy bond topologies for 
clusters of size N <1. Linear poly tetrahedral networks (b-e) and a 
closed 5-loop of face sharing tetrahedral (f) exhibit nonzero values of 
^ 2125 ^ 3235 ^ 434577545 , and/or 77555 . The octahedron (g) is the only 
maximally bonded cluster for A/ < 13 that is not composed of 
face-sharing tetrahedra. It has the same number of spheres and bonds 
as the 3-tetrahedron (d) and a different common neighbor signature, 
71200 = 2. 


favor the l^uid state to the point of removing the triple point 
from Fig. 0(a)l^. The square-well model also treats solvent 
in an implicit manner; explicitly accounting for solvent and the 
long-ranged hydrodynamics it mediates may b e imp ortant for 
colloidal crystallization under certain conditions^^lll^ 


those bonds (see Fig. (a)). We defined the relative num¬ 
ber Uabc = Nabc/N, where N is the number of particles. This 
analysis identifies gaseous, liquid, crystalline, and polytetrahe- 
dral configurations. Perfect face-centered cubic and hexago- 
nally close packed crystals exhibit nonzero values only of 77423 
(bulk hexagonal close-packed (hep)), 77424 (bulk hep and fee), 
77212, and 77312 (boundaries). Weakly interacting gases exhibit 
few bonds and uniformly low values of all common neighbor 
metrics, while weakly structured liquids exhibit large values of 
77200- Networks of face-sharing tetrahedra exhibit large values 
of 77323,77434,77545, and/or 77555 (see Fig.[^(b-f)). We character¬ 
ize the crystallinity by the fraction /c of particles participating 
in at least one 423 or 424 bond. 

To quantify the difference in crystallinity between dynamic 
simulations and simulations begun from fee or bcc crystals (see 
below), we use the ‘distance-to-equilibrium’ parameter 


Aeq = max 


/ _ ^random 

I ^42x ^42x 

I perfect 

V 


„bcc 

^666 


-n 


,random 
666 


perfect 

%66 


( 2 ) 


Here the superscript ‘random’ describes yield of common 
neighbor types from dynamic (randomly-initialized) simula¬ 
tions, the superscripts ‘fee’ and ‘bcc’ describe yields obtained 
from fee- and bcc-initialized simulations, respectively, and the 
superscript ‘perfect’ describes the yield of a given environ¬ 
ment in a bulk (perfect) fee or bcc crystal, i.e. = 6 and 


perfect _ . 

^666 — 

In Figs. [ 7 ] [8 [TO and jj_ we use the following color code 
to denote particle environments. Particles participating in crys¬ 
talline common neighbor environments (423 or 424) are col¬ 
ored green. Particles that do not, and that participate in poly- 
tetrahedral common neighbor environments (323, 434, 545, or 
555) are colored red. The remaining particles that participate 
in in liquid common neighbor environments (200) are colored 
blue. The remaining particles that participate in other common 
neighbor environments abc with a > 2 are colored magenta. 
The remaining (gas) particles are colored gray. 


1.2 Thermodynamics 


1 Methods 

1.1 Structure characterization 

We characterized the dynamics by performing a common 
neighbor analysis of the network of spheres linked by fa¬ 
vorable square-well interactions, similar to the analysis de¬ 
veloped by Honeycutt and Andersen®. At regular time in¬ 
tervals, we recorded the number of bonded pairs of parti¬ 
cles Nabc with a common neighbors, b bonds among those 
common neighbors, and c common neighbors participating in 


We tested the thermodynamic stability of finite-size crystals by 
performing virtual-move Monte Carlo simulations (see below) 
initialized from compact crystals containing approximately 
1000 spheres in a simulation box with an overall hard-core 
packing fraction of (j) =0.1. We used a 923-particle cuboc- 
tahedron for the fee crystal® and a 1001-particle cuboctahe- 
dron for the bcc crystal®. We initialized these crystals with 
nearest-neighbor distances J = 2(1 + A/2)/?o, placing nearest 
neighbors in the middle of their interaction range. We defined 
the boundary of finite-size fee stability (Fig.|^(a)) by the con¬ 
tour where 7742 ^ = 77423 + 77424 = 1 • (As shown by the relatively 
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Figure 3 (a) Translational and (b) rotational cluster diffusion constants, as a function of hydrodynamic radius R^, for virtual-move Monte 
Carlo simulations of tetrahedral clusters composed of infinitely attractive square-well spheres. We show results for various interaction ranges 
A (legend). The clusters used in these calculations were tetrahedral, and ranged in size from 1 to 120 particles; panel (c) shows a cluster of 120 
particles. Algorithm parameters were At = 4XRq, Aj- = 1, and pt = OA{RoAr/At)'^ (see text). Cluster diffusion constants approximate the 
Stokes solutions (A and A ^h~^)- For comparison, we show also the free-draining solutions (A and A which 

describe substantially slower collective motion. 


sharp decay of n 42 x in Fig.[^(a), the location of the boundary is 
relatively insensitive to choice of threshold.) 


We calculated the boundary of stability of bulk fee crystals 
by performing single-particle Monte Carlo (SPMC) direct co¬ 
existence simulations of 1000 spheres in a slab geometry at 
various temperatures and interaction ranges, choosing over¬ 
all packing fractions that allowed sufficient sampling of both 
fluid and crystal phases. We initialized the crystal slabs with 
nearest-neighbor distances J = 2(1 + A/2)/?o, and we initial¬ 
ized the fluid phases with random conflgurations without hard¬ 
core overlaps. As shown for example by the black points in 
Fig. these simulations allowed us to determine the coexis¬ 
tence concentrations for the fluid phase (gas or liquid, depend¬ 
ing on T and A) and the crystal phase (not shown). We deflned 
the boundary of bulk crystal stability at 0 =0.1 for each A as 
the temperature at where the interpolated fluid coexistence con¬ 
centration (black curves in Fig.[^ intersects 0 = 0.1. 


We calculated the boundary of stability of the bulk liquid by 
performing SPMC Gibbs ensemble Monte Carlo simulations of 
1000 spheres separated in two boxes that exchange spheres and 
volumeEH. Analogous to the direct coexistence simulations, we 
performed the Gibbs ensemble simulations at a range of tem¬ 
peratures for each interaction range, with overall packing frac¬ 
tions chosen to allow sufficient sampling of both phases, and 
we determined the boundary of liquid (meta)stability at 0 = 0.1 
by interpolating the fluid coexistence curves (see blue points 
and curves in Fig. [^. We initialized both the gas and liquid 
box with random configurations without hard-core overlaps. 


1.3 Dynamics 


To approximate the overdamped dynamics of strongly- 
associating particles in solution we used the virtual-move 
Monte Carlo algorithm®! (specifically, the version of the al¬ 
gorithm described in the appendix of Ref.^D). Under this algo¬ 
rithm, which satisfies detailed balance, particles move locally 
according to the gradients of potential energy they experience, 
and collectively with a rate that can be controlled to a degree by 
the user. We parameterized the algorithm in a manner similar 
to that described in Ref.®!, in order to ensure that tightly-bound 
clusters of particles of hydrodynamic radius Rr diffused with 
rates close to those predicted by the Stokes’ law, 


A = 
A = 


6Kr{R^ ’ 
k^T 

8;rT7RH^ 


(3) 


The natural time unit of this motion is then the Brownian time 
scale 


t](2/?o)2 

^0 = j rj. 

k^T 


(4) 


where r] is the (implicit) solvent viscosity and Ab T is the ther¬ 
mal energy. Simulations performed at particular values of 
k^Tje and A can therefore be considered to apply to a wide 
range of absolute particle sizes Rq, the latter determining only 
the value of ^o- For example, for spherical colloidal particles 
with radius Rq = 50 nm at room temperature {T = 293 K) 
in water (r] = 1.00 x 10“^ Pa s), the Brownian timescale is 
^0 = 2.5 X 10-4 s. 

In the Appendix we describe in detail the procedure we used. 
Briefly, the virtual-move algorithm generates collective Monte 
Carlo moves by proposing trial ‘virtual’ particle translations or 
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2 RESULTS 


rotations, and probabilistically recruiting neighboring particles 
to join this motion, in an iterative fashion. The resulting trial 
move is accepted with a probability ensuring detailed balance. 
In addition, one is free to attenuate the rate at which collec¬ 
tive motion is accepted, by imposing what are effectively ki¬ 
netic constraints. We chose these constraints in order to enforce 
Eq. 0. 

Each Monte Carlo move begins with either a trial translation 
or a trial rotation, chosen with probability pi and = I — Pu 
respectively. For translations, we randomly selected a particle 
and translated it randomly within a ball of radius At. For ro¬ 
tations, we randomly selected a particle, randomly selected a 
second particle within the interaction range of the first, and ro¬ 
tated the second particle by an angle, chosen uniformly from 
the range (—Ar, Ar), around a randomly-oriented axis h passing 
through the center of the first particle. We chose At = 4A/?o, 
Ar = 1, and pi = 0.4(7?oAr/At)^. As discussed in the Ap¬ 
pendix, we found that with this choice of parameters we could 
enforce Eq. 0 by suppressing the acceptance rate for transla¬ 
tion and rotation of a cluster of N particles of hydrodynamic 
radius 7 ?h by factors N~^R^ and respectively. We 

chose to maximize the ratio of rates of internal cluster relax¬ 
ation to whole-cluster diffusion, by working with the smallest 
trial displacement At that is large enough to induce substantial 
collective motion (i.e. is large enough to ensure that Stokes’ 
law could be maintained). This choice is somewhat arbitrary, 
and whether it is physically appropriate will likely depend on 
details of the experimental system one wishes to model, but we 
note that one has some freedom to infiuence this ratio if neces¬ 
sary. 

In Fig. we show measured diffusion constants for tetra¬ 
hedral clusters of between 1 and 120 square-well spheres, in 
the ^ 0 limit. and approximate the Stokes so¬ 

lutions (Eq. over a broad range of cluster sizes and interac¬ 
tion ranges. In this respect our procedure therefore captures 
an important aspect of solvent-mediated diffusion, without rep¬ 
resenting solvent explicitly. Note that ‘long-ranged’ hydrody¬ 
namic coupling03 is not taken into account by this pr ocedu re; 
to do so, one should represent solvent more explicitlyE^ESEl, 
Simple implementations of Brownian (Langevin) dynamics in¬ 
tegrators, and single-particle Monte Carlo simulations in the 
limit of zero trial displacementlSl, result instead in the ‘free- 
draining’ behavior cx and oc R^U^. As shown in 

Fig. (note the logarithmic scale), such diffusion is signifi¬ 
cantly slower than Stokes’ diffusion, even for relatively modest 
cluster sizes. 

As discussed in the Appendix, our procedure yields a time 
per Monte Carlo cycle of 

6 .2 

^cycle — (5) 

where to is the Brownian time scale (Eq.|^. We present results 


relative to the physical time unit ^o- 

We carried out simulations of 1000 square-well spheres, in 
periodically-replicated cubic simulation boxes, at a hard-core 
packing fraction of (j) =0.1. We carried out independent sim¬ 
ulations for interaction ranges between (and including) the val¬ 
ues A = 0.005 and 1.35, and for temperatures ranging from 
kBT/£ = 0.06 to 0.86. We initialized dynamic simulations with 
random configurations, under the constraint that the particle 
hard cores could not overlap (equivalent to equilibrium con¬ 
figurations m\hQ k^T/£ ^ oo limit). 

An open-source library for implementing the virtual- 
move Monte Carlo algorithm is available at http://vmmc.xyzPSl. 

2 Results 

2.1 Dynamic and thermodynamic phase dia¬ 
grams 

In Fig. Ua) we show in the temperature-range plane the fee 
crystal yield /c seen in dynamic simulations at time t = 10^to’, 
panel (b) shows yield also at times lO^^o and lO^^o- High yield 
(green) is found in a localized region of parameter space. Note 
that the phase diagrams of Fig. inter sect the diagram of Fig. 
via three vertical lines corresponding to particular values A = 
0.18, A =0.25, and A =0.33. 

A necessary condition for high crystal yield is that the fee 
crystal is stable thermodynamically; this condition holds for 
our compact fee crystals below the solid green curve, marked 
“g-c (finite)” in Fig. Note that the dashed green bulk gas- 
crystal coexistence curve (labeled “g-c”) derived from direct 
coexistence simulations in a slab geometry lies above the finite- 
size curve. This difference simply refiects the fact that a fi¬ 
nite crystalline cluster with free boundaries can melt within 
the regime of bulk crystal stability, i.e. can be smaller than 
the critical cluster size. Note also that the solid green curve 
bends toward small A and small k^T at around A =0.32, 
k^T{£ = 0.44. This bend occurs because the fee crystal be¬ 
comes unstable with respect to to a body-centered cubic (bcc) 
crystal at large A and small k^T{£ (see below). 

A second necessary condition for high crystal yield expected 
from previous workESGSl is that a state point must lie within 
the regime of metastable liquid-gas phase coexistence. For 
large interaction ranges A we determined the boundary of liq¬ 
uid stability in bulk from Gibbs ensemble simulations (dashed 
gray curve in Fig. |^, a technique that eliminates interfacial 
effects by putting gas and liquid phases in separate boxes 
that interchange both volume and particlesEH. For interaction 
ranges smaller than A =0.15, we could not accurately deter¬ 
mine gas-liquid coexistence because crystallization rapidly oc¬ 
curred within the liquid box. Instead, we estimated the onset 
of transient liquid-like structure as the curve below which our 
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2 .1 Dynamic and thermodynamic phase diagrams 


2 RESULTS 



Figure 4 (a) Crystal yield fc from dynamic simulations after t = 10^to as a function of interaction range A and temperature k^T/e, for 
systems of 1000 square-well spheres at hard-core packing fraction 0 = 0.1. The dashed (solid) green curve indicates the boundary of stability 
of bulk (finite-size) fee crystals, and the dashed (solid) gray curve indicates the boundary of stability of the bulk (finite-size) liquid (see text 
for details). Symbols indicate representative state points for the various dynamic regimes shown in Figs.[^|^[^[^ and[TT] (b) Yield at three 
times. 
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2 RESULTS 


2.1 Dynamic and thermodynamic phase diagrams 



Figure 5 (a) Common neighbor metric n 42 x indicating fee crystallinity after simulations of length t = 10^to initiated from perfect fee crystals, 
(b) Common neighbor metric ^666 indicating bee crystallinity after simulations of length t = 10^ to initiated from perfect bee crystals, (c) 
Distance-to-equilibrium Agq (defined in Section p~T| indicating how close the common neighbor metrics in the dynamic simulations (Fig.|^ 
get to the closest of the two crystal-initiated simulations. 


dynamic simulations attained a relative number of liquid-like 
bonds tZ2oo > 0.1 at some point during our dynamic simulations 
(solid gray curve in Fig. see also Fig. 12 (g)). We find that 
this curve coincides with the bulk liquid curve over the interval 
of interaction ranges for which both could be calculated, indi¬ 
cating that metastability of the liquid is not strongly influenced 
by the existence of free boundaries for system sizes on the or¬ 
der of 1000 particles. Although it is not clear to what extent the 
liquids at small interaction range can be considered metastable, 
comparison of the onset of high crystal yield (green pixels) with 
the liquid boundary (dashed and solid curves) is consistent with 
crystallization coinciding with the onset of transient liquid or¬ 
der and/or an extrapolation of the metastable liquid curve. 


While our results show that the stability of the crystal and 
the onset of transient liquid structure are necessary conditions 
for rapid crystallization, they are clearly not sufficient: large 
regions of parameter space below the crystal and liquid bound¬ 
aries appear blue or red in Fig. indicating low crystal yield 
after t = 10^to, despite the fact that the systems must eventually 
assemble into a thermodynamically favored fee crystal. The 
eventual (infinite-time) fate of the system can be inferred from 
Fig. in which we plot in panels (a) and (b) the crystal yield 
that results from simulations initiated from a single fee or bcc 
crystal, respectively. Note that the fee crystal becomes unstable 
with respect to the body-centered cubic (bcc) crystal at large A 
and low k^Tje, because bcc spheres can accommodate 8 near¬ 
est neighbors and 8 second-nearest neighbors at these values of 
A, while fee spheres can only accommodate 12 nearest neigh¬ 
bors. In this region of parameter space the bcc crystal is there¬ 
fore lower in energy than the fee crystal. Panel (c) displays a 
‘distance-to-equilibrium’ parameter (see Section pTT] ) that sum- 


■ Metastable gas 



t/to 


Figure 6 (a) Snapshot from a simulation with A = 0.03 and 
k^Tje = 0.32. The system remains in the metastable gas phas e up to 
time t = Particle color code is described in Section |l.l| (d) 

Time series of the crystal yield /c (green), liquid (blue) and 
poly tetrahedral common neighbor metrics. 


marizes how close dynamic simulations come to equilibrium: 
anything not shown green corresponds either to a metastable 
liquid or to an arrested gel. 

As we will discuss below, low crystal yield within the ther¬ 
modynamically stable crystal region is due to one of two kinetic 
effects, depending on the state point: either nucleation from the 
liquid is slow, or crystallization is arrested by gelation. Further¬ 
more, we find that the kinetics of crystallization varies strongly 
even within the region of high yield: crystallization may pro¬ 
ceed in a two-step pathway via a liquid-like intermediate, or it 
may proceed directly from a relatively homogeneous gas. The 
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2.4 Two-step crystallization 
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• One-step 



Figure 7 (a-c) Snapshots from a simulation with A = 0.03 and k^T/s = 03 that exhibits one-step crystallization (a) before nucleation, (b) 
after nucleation, and (c) at the end of the simulation (t = 10^ to)- (d) Time series of the crystal yield /c and liquid and poly tetrahedral common 
neighbor metrics. Liquid-like environments (blue curve and particles) are seen throughout crystallization, but the crystal that nucleates and 
grows does not have substantial liquid-like character. Arrows indicate the time points of the snapshots. 


following five subsections discuss the five qualitatively differ¬ 
ent dynamic regimes encountered when broadly varying the in¬ 
teraction range and strength. 

2.2 Metastable gas 

The black square in Fig. lies in the metastable gas regime. 
Here, as shown in Fig.[^ the system remains in a gas state, with 
low values of all common neighbor metrics. In this regime the 
face-centered cubic (fee) crystal is the stable state, while the 
liquid is unstable with respect to the gas. The metastability of 
the gas for times up to ^ = 10^ indicates the existence of large 
free-energy barriers for direct crystal nucleation from the gas. 

2.3 One-step crystallization 

As the temperature decreases below the metastable liquid tran¬ 
sition (dashed and solid gray curves in Fig. the crystal yield 
increases, as indicated by the sharp change from red (low yield) 
to green (high yield) pixels in Fig. As shown in Fig. 
(black circle in Fig. |^, dynamic crystallization pathways near 
the metastable liquid transition involve crystal nucleation from 
a fluid with substantial liquid-like fiuctuations but no signifi¬ 
cant gas-liquid phase separation. Before the sharp increase in 
crystal fraction fc at around 6 x 10^ (Fig- 0(b)), there are sub¬ 
stantial fiuctuations in the liquid common neighbor metric ^ 200 , 
but the average value of ^200 remains less than 0.2 (one liquid¬ 
like bond per five spheres). The nucleation event at around 
6 X lO^^o does not occur at the expense of liquid-like structure, 
as it would if the crystal nucleated from within a liquid-like 
droplet. Instead, the value of ^200 increases during the nucle¬ 
ation event. Thus, Fig. [^illustrates a one-step pathway that ap¬ 
pears to be facilitated by strong but non-critical density fiuctua¬ 


tions. A similar pathway was identified in RefH As discussed 
below, the region of parameter space in which we observe a 
one-step pathway is very narrow, consistent with the fact that it 
was not found in many other studies. 

2.4 Two-step crystallization 

Beyond this narrow region of one-step nucleation we find a 
broad range of parameters where crystallization occurs via a 
two-step pathway, illustrated in FigJ^(a-d) and (e-g) (black star 
and diamond, respectively, in Fig. [^ First, liquid-like clusters 
(Fig.[^(a) and (e)) quickly nucleate, grow, and merge, result¬ 
ing in an increase in the liquid common neighbor metric ^ 200 - 
Later, crystals (Fig. (b) and (f)) nucleate from within those 
droplets, resulting in a decrease in ^200 and an increase in the 
crystallinity metric /c. The nucleation time increases with in¬ 
creasing range A, as can be seen by comparing the common 
neighbor time series of Fig. [^(d) and (g). For A > 0.2 the time 
for crystal nucleation from the liquid exceeds our simulation 
time t = lO^^o- 

To illustrate the crossover from one-step to two-step path¬ 
way more generally, we show in Fig. [^a) a parametric plot 
of the liquid common neighbor metric ^200 versus the frac¬ 
tional crystal yield /c, for a slice of state points (A = 0.03 and 
0.22 <k^T/£< 0.3) having crystal yields /c > 0.7 ait = lO^^o- 
Most of these systems {k^T/e < 0.28, including the example 
k^T{£ = 0.28 from Fig.[^ follow a pronounced two-step path¬ 
way. First, ^200 increases with little increase in /c, correspond¬ 
ing to liquid droplet nucleation, growth, and coalescence. Sub¬ 
sequently, ^200 decreases and /c increases, corresponding to 
crystal nucleation (fast or slow) from within the liquid droplet. 
In contrast, the example system with AbT/e = 0.3 (Fig.j^ does 
not exhibit pronounced two-step nucleation; instead, a crystal 
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2 RESULTS 


2.4 Two-step crystallization 



Figure 8 (a-c) Snapshots from a simulation with A = 0.03 and k^T /e = 0.28 exhibiting fast two-step crystallization: (a) before crystal 
nucleation, (b) after crystal nucleation, and (c) at the end of the simulation {t = 10^to)- (d) Time series of the crystal yield /c and liquid and 
polytetrahedral common neighbor metrics. Arrows indicate the time points of the snapshots, (e-f) Snapshots from a simulation with A = 0.17 
and k^T/£ = 0.48 exhibiting slow two-step crystallization: (e) before crystal nucleation and (f) after crystal nucleation, at the end of the 
simulation {t = lO^to)- (g) Time series of the crystal yield /c and common neighbor metrics. 
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2 RESULTS 



^200 A 


Figure 9 Parametric pathway diagram illustrating the evolution of 
liquid structure (^ 200 ) on the horizontal axis and crystalline structure 
(/c) on the vertical axis, for a slice of state points with A = 0.03 and 
0.22 <kBT/£ < 0.3. All trajectories show high crystal yield /c > 0.7 
after t = 10^ to- For most state points (k^T/s < 0.28) crystallization 
follows a two-step pathway, where first ^200 increases and then ^200 
decreases while /c increases. For k^T/s = 0.3 crystallization 
proceeds via largely a one-step mechanism, with n 2 oo remaining low 
throughout assembly, (b) Maximum value of ^200 during assembly as 
a function of A and k^T/e, restricted to those state points for which 
/c > 0.7 after t = 10^ to- Most points are green, indicating a two-step 
pathway with large intermediate values of ^ 200 - A narrow strip of 
state points near the extrapolated location of the gas-liquid curve are 
blue, indicating one-step assembly for which n 2 oo is low throughout. 


nucleates directly from the gas. 

In Fig.[^(b) we show the maximum value of ^200 along each 
pathway for which /c > 0.7 after t = lO^^o- The narrow strip of 
values along the top of this region, near the metastable liquid 
transition, exhibit one-step behavior, with correspondingly low 
(blue) values of ^200 . 


2.5 Metastable liquid 


At larger interaction ranges (A > 0.2) and for temperatures 
close to the gas-liquid curve we find that the liquid remains 
metastable up to times t = 10 ^to. The liquid metric ^200 is large, 
and the crystal metric /c increases until it reaches a plateau that 
persists until the end of the simulation. An example trajectory 
is shown in Fig. 10 (black pentagon in Fig.[^. For larger sys¬ 
tems and times that are longer (but still accessible to the cor¬ 
responding experiments), this region of parameter space may 
give rise to good crystals. 


2.6 Gelation 

At low temperature, our simulated systems dis play t he fast 
formation and persistence of polytetrahedral gels F^I^^ I, which 
are bonded networks of face-sharing tetrahedra. As illus¬ 
trated for small networks in Fig. (b-f), polytetrahedral net¬ 
works exhibit nonzero values of the common neighbor metrics 


# Metastable liquid 



Figure 10 (a) Snapshot from a system with A = 0.2 and 
k^Tje = 0.54 which remains as a metastable liquid up to time 
t = 10^to- (b) Time series of the crystal yield /c and liquid and 
polytetrahedral common neighbor metrics. 


t^2i25t2323,t2434,t2545, and ^2555. All but the 27212 metric do not 
appear in perfect close-packed crystals; 212 environments ap¬ 
pear on the 100 and 110 surfaces of fee crystals. As shown 
for example in Fig. 11 (black triangle in Fig. the metrics 
^323,^7434,77545, and/or 77555 increase quickly and remain large 
up to times t = lO^to. while metrics characterizing crystals (/c) 
and mobile liquids (77200) remain low. 

We find that the decrease in crystal yield at low tempera¬ 
ture is accompanied by proliferation of these polytetrahedral 
networks. This correspondence can be seen across parameter 
space by referring to Fig.[^(a-e), in which we compare (a) the 
crystal yield fc after t = l^o with (b-e) the maximum value of 
the polytetrahedral common neighbor metrics 77323,77434,77545, 
and 77555 achieved during the simulations. Although each com¬ 
mon neighbor metric shows somewhat different dependence on 
A and AeT/e, comparison of Fig. 12 (a) with Fig. 12 (b-e) 
shows that pathways involving large maximum values of the 
polytetrahedral common neighbor metrics account for most of 
the region in Fig. 12 (a) where the close-packed crystals are 
stable but yield is low. 


The predominance of polytetrahedral common neighbor 
metrics at low temperature suggests an explanation for the dy¬ 
namic inaccessibility of the bcc crystal. Since the region of 
parameter space where the bcc crystal is stable (Fig. (b)) is 
contained within the region where the polytetrahedral common 
neighbor metrics are large (Fig.[^(a-e)), there is no region of 
parameter space where bcc crystallization proceeds efficiently. 
Rather, polytetrahedral gelation forestalls bcc crystallization 
wherever the bcc crystal is thermodynamically stable. 


The only large region below the gas-liquid and liquid-crystal 
curves where neither /c nor the polytetrahedral metrics are 
large is the region marked “metastable liquid” in Fig.[T^(a). As 
discussed above, crystal yield is low in this regime because the 
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2 RESULTS 


2.6 Gelation 




Figure 11 Snapshots from a simulation with A = 0.27 and k^T/e = 03 that forms a polytetrahedral gel (a) before gelation, (b) after gelation, 
and (c) at the end of the simulation (t = 10^to)- (d) Time series of the crystal yield fc and liquid and polytetrahedral common neighbor 
metrics. Arrows indicate the time points of the snapshots. 



Figure 12 (a) Relative crystal yield fc after t = lO^tQ as a function of interaction range A and temperature k^T/s. (b-e) Maximum values of 
the polytetrahedral common neighbor metrics (b) W 323 , (c) W 434 , (d) W 545 , and (e) W 555 over the course of the simulations, as a function of 
interaction range A and temperature k^Tje. The region of parameter space in the lower right of (a) (“gel”) where fc <0.1 is mostly 
associated with large values of ^ 323 ^,^ 43 ^ 4 ’^^54?’ and/or (b-e), except for the region labeled “metastable liquid”, (f) Nucleation time tnuc 
(time to first achieve fc = 0.7) on a logarithmic scale, (g) Maximum values of the liquid common neighbor metric ^ 200 - 
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nucleation time is longer than not because there is sub¬ 

stantial polytetrahedral gelation. To see that nucleation times 
extrapolate beyond 10^in the metastable liquid regime, in 
Fig.[T^(f) we plot the time (on a log scale) at which each sys¬ 
tem first achieves a crystal yield of fc = 0.7. Nucleation times 
increase to the right and top within the high-yield region until 
they pass beyond the t = 10^to window. Fig. 12 (g) shows that 
beyond the high-yield region, the maximum value of ^200 re¬ 
mains high, indicating that these systems achieve similar levels 
of liquid structure as in the high-yield region, the main differ¬ 
ence being that crystals have not yet nucleated from the liquid 
after lO^^o- 


3 Conclusions 


Our results confirm that colloidal crystallization displays fea¬ 
tures common to many examples of self-assembly, in that it 
happens effici ently in only a small regime or ‘sweet spot’ of pa¬ 
rameter spacelSElHISl. We confirm that efficient crystal nucle¬ 
ation of spherical particles with short-range attractions happens 
via a two-step pathway throughout most of parameter space. 
However, by performing a systematic investigation of colloidal 
crystallization as a function of interaction strength and range, 
we have found that largely one-step crystallization from the gas 
can occur for certain combinations of parameters. Specifically, 
over a narrow range of interaction strengths and ranges near 
the extrapolated location of the metastable liquid-gas boundary, 
we find that crystallization occurs without significant accumu¬ 
lation of liquid-like order within the relatively short time frame 
oft = 10^r]Rl/kBT. 

Our results show that low crystal yield at low temperatures 
or large interaction strengths is accompanied by particular local 
bond geometries. Crystal yield is low whenever there are many 
common neighbor configurations associated with face-sharing 
tetrahedra. Royall et al. showed that gelation in a colloid- 
depletant mixture can be explained by formation of overlap¬ 
ping networks of locally-favored state s, eac h defined as a 
maximally-bonded cluster of size m <13111601 

For short-range 

interactions as explored in Ref.®l, most locally-favored states 
are sections of a 13-particle icosahedron consisting of 20 face¬ 
sharing tetrahedra, and are therefore clusters of face-sharing 
tetrahedra®. The only exception is the octahedron (Fig.|^(g)), 
which has the same number of spheres (6) and bonds (12) as 
three face-sharing tetrahedra (Fig.[^(d)), but is not composed 
of tetrahedra. The octahedron thus shows a distinct common 
neighbor signature: instead of exhibiting nonzero values of the 
polytetrahedral common neighbor metrics, the octahedron ex¬ 
hibits only a nonzero value of the common neighbor metric ^200 


that is prevalent at temperatures above gelation (Fig. (g)). 
Our findings are largely consistent with the results of Royall 
et al ., suggesting that key features of nonequilibrium colloidal 


assemblies can be captured by the square-well model. 

We note also that within the square-well system, gelation can 
be detected by the study of local bond environments and does 
require identifying maximally-bonded clusters. Indeed, our re¬ 
sults suggest that there may be some local configurations seen 
in gels that do not participate in maximally bonded clusters: 
while 323, 434, and 555 bonds associated with polytetrahedral 
gelation are found in maximally bonded clusters, 545 bonds are 
not. Instead, 545 bonds are found in curved, linear polytetra¬ 
hedral motifs as shown in Fig. (e). As seen by comparing 
Fig.|^(e) to Fig.[^(f), these motifs do not maximize the num¬ 
ber of bonds because they do not close into complete loops of 
face-sharing tetrahedra. 

Finally, our dynamic protocol suggests which combina¬ 
tions of temperature and colloid interaction range will yield 
best crystallization on the time scale lO^^o, where = 
7](2Ro)^/I cbT. For example, for spherical colloidal particles 
with radius Rq = 50 nm, at room temperature (T = 293 K) in 
water (r/ = 1.00 x 10“^ Pa s), our results predict that crystal¬ 
lization after r = lO^^o = 25 s will be best for an attractive in¬ 
teraction of range 0.04Ro = 2 nm and strength O.SIIcbT = 0.19 
kcal/mol. 
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5 Appendix: Virtual-move Monte 

Carlo algorithm 


To efficiently approximate the overdamped and hydrodynami- 
cally coupled dynamics of strongly associating particles in so- 
lutio n, we used the virtual-move Monte Carlo (VMMC) algo¬ 
rithm®®, specifically the variant described in the appendix of 
Ref.®. We parameterized the algorithm to satisfy the Stokes 
solutions for the translational and rotational diffusion of clus¬ 
ters of hydrodynamic radius Ru, 


A = 

Dr = 


kBT 

6;rT7RH ’ 
kBT 

SttviRu^ ’ 


( 6 ) 


while allowing as much internal relaxation of each cluster as 
possible. Our parameterization allows us to present results 
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5 .1 Review of the VMMC algorithm 


across interaction ranges and particles sizes relative to the nat¬ 
ural Brownian time unit 

to = n{2Ro?/kBT, (7) 

where rj is the solvent viscosity, Rq is the hard-core radius of 
the spherical particles, and k^T is the thermal energy. 

The VMMC algorithmE3E2 moves individual particles and 
groups of particles with attempt and success frequencies de¬ 
signed to (1) preserve the correct equilibrium distribution, (2) 
ensure that particles move according to gradients in the po¬ 
tential energy, and (3) allow the dependence of diffusion co¬ 
efficients on cluster size (e.g. Eq. (|^) to be controlled. 
The algorithm achieves this by proposing individual Monte 
Carlo moves, self-consistently generating individual or col¬ 
lective moves from the proposed moves, and accepting those 
moves in a way that satisfies the above three conditions. 

In our implementation, trial individual translations are at¬ 
tempted with probability pi by randomly selecting a particle 
and then attempting to translate it randomly within a ball of 
radius At. Because our particles are spherically symmetric, col¬ 
lective rotations cannot be generated from trial rotations about 
the center of mass of a single particle. Our trial rotations 
are attempted with probability pj^ = 1 — pt by randomly se¬ 
lecting a particle, randomly selecting a second particle within 
the interaction range of the first, and then attempting to rotate 
the second particles by a randomly-chosen angle in the range 
(—Ar,Ar) around a randomly oriented axis h centered at the 
first particle. 

Acceptance rates in the VMMC algorithm consist of three 
factors: (1) a factor built on the Metropolis criterion ensuring 
that the system relaxes toward equilibrium and particles move 
according to gradients in the potential energy (2) a factor ensur¬ 
ing that motions of clusters are not oversampled with respect to 
the motion of isolated particles, and (3) a factor enforcing a 
prescribed dependence of diffusion coefficients on cluster size. 
The first factor is generic for any application of the VMMC 
algorithm. The second factor is usually used to produce real¬ 
istic dynamics, but can be omitted when the algorithm is used 
only to sample an equilibrium distribution. The third factor has 
a general form that depends on the scaling of diffusion coeffi¬ 
cients D[ and Dj- on hydrodynamic radius /?h (e.g. A ^ 
and A oc in Eq. [^, but the parameters of the algorithm 
{pi. At, and Ar) must be tuned to ensure that the prefactor of 
the diffusion laws match the prescribed values® Below, we 
review the VMMC algorithm and discuss our parameter opti¬ 
mization. 

5.1 Review of the VMMC algorithm 

We followed the ‘symm etrize d’ version of the VMMC algo¬ 
rithm discussed in Refs.®®. Trial ‘virtual’ single-particle 


translations and two-particle rotations are generated as dis¬ 
cussed above, with probability pi and pr = pt —I, respectively. 
Trial collective moves are generated by iteratively testing each 
link between particles inside the moving group and particles 
outside the moving group, starting with an initial translation or 
rotation, until no links remain to be tested. When each link i j 
is tested, the particle outside the moving group (j) is pre-linked 
to the moving group with probability 

= ^j-max (0,1 - exp (j 3 (t/^ - Ui^j ))), (8) 

where = 1 if particles i and j are within their mutual inter¬ 
action range and = 0 otherwise, Uij is the initial interaction 
energy, and A'y is the interaction energy when particle i exe¬ 
cutes a virtual move but particle j does not. After making a 
virtual move particle i is returned to its initial coordinates. If 
particle j is pre-linked to the moving group, then the probabil¬ 
ity of the reverse move (translation or rotation in the opposite 
direction) is calculated, 

^reverse ^ j _ ^ ^ 

where interaction energy when particle i does the 

reverse virtual move and particle j stays still. After the reverse 
virtual move i is restored to its initial position. Then, the pre¬ 
linked particle j is linked (admitted) to the moving group with 
probability 

Otherwise, the link is marked as frustrated. This procedure is 
continued until it converges on a final moving group, with all 
links having been tested between particles inside and outside 
the group. 

Since the algorithm will move the entire group by the ini¬ 
tial translation or rotation, the acceptance of collective moves 
must be scaled by a factor of 1/A to prevent over-sampling the 
motion of particles that are in large clusters. This is achieved 
by rejecting moves in situ when the number of particles in a 
moving group grows larger than \/x (translations) or 2/x (rota¬ 
tions), where v is a random number on the interval (0,1] chosen 
at the beginning of each Monte Carlo step. The factor 2 appears 
because only rotations of groups of at least two particles are ex¬ 
plicitly simulated for spherical particles; rotations of individual 
particles can be assumed to occur at any arbitrary rate while 
having no effect on the center of mass motion of the particles. 

An additional scaling of acceptance probabilities can be ap¬ 
plied to control the dependence of diffusion on cluster size. The 
Stokes scalings (Eq. can be perfectly enforced in the limit 
where particles are tightly bound to each other within well- 
defined clusters and do not interact outside of these clusters. 
In this limit, any trial Monte Carlo move results in the recruit¬ 
ment of the entire cluster into a moving group, followed by the 
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translation or rotation of the entire cluster. The VMMC algo¬ 
rithm can enforce Eq.j^in this limit by rejecting moves in situ 
when the hydrodynamic radius Ru of the moving group exceeds 
Rmin/y (translations) or Rmin/y^ (rotations), where R^i^ is the 
minimum possible hydrodynamic radius (see below) and y is 
another random number on the interval (0,1]. Following pre¬ 
vious implementations of the VMMC algorithmSSl we estimate 
the hydrodynamic radius Ru of a group ^ as a generalization 
of the radius of gyration, 

Rll = 10(1 (r — Tcenter) ^1 (H) 


where reenter is the group’s center of mass (center of rotation) 
and h is the direction of the translation (axis of rotation) for 
translations (rotations). This factor is the same for for trans¬ 
lations that occur in opposite directions, and for rotations that 
occur with opposite sense, as is required for detailed balance. 
We take r G ^ to include all points within the hard cores of the 
particles. The minimum hydrodynamic radius for both transla¬ 
tions and rotations (single-sphere translations or effective two- 
sphere rotations about an axis h parallel to the separation vector 
between the particles) is the physical sphere radius Rq (it would 
be Rq / without the factor of 10 that appears in Eq. (11)). 

Once a moving group has been generated that is not rejected 
in situ due to its number of particles or hydrodynamic radius, 
two additional factors contribute to its acceptance probability. 
First, the move is rejected if there are any frustrated links be¬ 
tween particles inside and outside the moving group. Second, 
moves that remain valid are accepted with probability 


Wacc = min I 1 , jq exp (-j3 - Uij)) 1 , (12) 

\ (O)oop / 

where the product runs over all pairs of particles (/ in the mov¬ 
ing group and j outside it) that are non-interacting before the 
move and have positive pair energy after the move, or vice 
versa. Together, these factors ensure that the system satisfies 
superdetailed balance, a condition that implies detailed bal- 
anceED. For square well spheres that have only zero, negative, 
or infinite positive interaction, Eq. [^reduces to a rejection if 
the move results in any hard-core overlap, 

Wacc= n 0 ('•,•' 7 - 2 ^ 0 ), (13) 

where 6 is the heaviside step function. 

For rotations of moving groups large enough to interact with 
their periodic images we imposed an additional rejection if a 
move resulted in a hard core overlap with a periodic image. 
Such overlaps occurred only for gels. 


5.2 Parameter optimization 

The prescribed dependence of diffusion coefficients with clus¬ 
ter size (Eq. is derived for the VMMC algorithm in the limit 


of vanishingly narrow potential energy wells. In this limit, trial 
moves always take particles out these wells. If the wells are 
deep relative to ke T and the clusters are isolated, this causes the 
algorithm to recruit the entire cluster into the moving group, al¬ 
ways accept the move, and thus generate diffusion coefficients 
dictated by the in situ rejection of collective moves as a func¬ 
tion of hydrodynamic radius. 

In practice, potential energy wells are not vanishingly nar¬ 
row. Properly modeling the motion of such systems requires 
(1) allowing degrees of freedom internal to the clusters to relax 
and (2) ensuring that the prescribed diffusion laws are obeyed 
even when whole-cluster moves are not always generated. We 
sought to satisfy these conditions by choosing algorithm pa¬ 
rameters that allow clusters to internally relax as much as pos¬ 
sible without violating the Stokes scaling. We achieved this by 
selecting the smallest values of At and Ar that resulted in the 
Stokes solutions for a test set of tetrahedral clusters of size 1 
to 120 (8 spheres along each edge) for a range of interaction 
ranges A and T ^0 (infinite square well attractive interaction). 
We determined the optimal algorithm parameters in four steps. 

First, we recorded the translational diffusion coefficient Dt 
from translation-only VMMC simulations of isolated tetrahe¬ 
dral clusters with range A =0.11 and compared it to the pre¬ 
dicted value from Eq. 


nStokes =D° — , 

‘ 'Rh 


where 




lOt 


cycle 


(14) 


(15) 


is the translational diffusion coefficient for a single sphere 
translating with Monte Carlo dynamics. As shown in Fig. 

(a), we found that for a range of tetrahedral cluster sizes 
A « for At < 0.2Ro, A - for A, = 0.4/?o, 

and A — for At > Ro- Although the translational dif- 

fusion coefficient is not quite saturated at the Stokes limit at 
At = 0.4/?o, we chose to parameterize At near this value be¬ 
cause this balances a diffusion coefficient near the Stokes limit 
with ample internal relaxation: as shown in Fig.[^(b), more 
single-particle than whole-tetrahedron moves are accepted for 
At = 0.4/?o, while the reverse is true for Ai= Rq. 

Second, we recorded the rotational diffusion coefficient Z)r 
from rotation-only VMMC simulations of the same tetrahedral 
clusters and compared it to the predicted value from Eq. 


^Stokes ^ ^0 



where 




ems(Ar)Ar^ 

6/cycle 


(16) 


(17) 
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Rh/Ro Rh/Ro 


Figure 13 (a) Ratio of translational diffusion coefficient to the value 
predicted by the Stokes solution (Eq. [EJ for translation-only VMMC 
simulations of tetrahedral clusters composed of spheres with infinite 
square well attractive interactions of range A = 0.11. (b) Ratio of 
accepted single-sphere moves to whole-cluster moves for the same 
set of simulations, (c) Ratio of rotational diffusion coefficient to the 
value predicted by the Stokes solution (Eq. for rotation-only 
VMMC simulations of the same tetrahedral clusters, (d) Ratio of 
translational to rotational diffusion coefficients for a four-sphere 
tetrahedron with infinite square well attractive interactions as a 
function of the attempt probability for translations. The ratio is 
normalized such that the Stokes solution corresponds to a value of 1 
(see Eq.|^. (e) Translational and (f) rotational diffusion coefficients 
vs hydrodynamic radius for simulations of tetrahedral clusters of 
various sizes composed of infinite square well attractive spheres with 
various interaction ranges A (legend), using the full translating and 
rotating VMMC algorithm with parameters At = 4ARo^ Ar = 1, and 
Pt = 0.4(Ro^r/At)^. Eor comparison, the Stokes solutions 
(Dt oc and oc and free draining solutions (A oc Rii~^ 

and Dr ^ Ru~^) are shown as lines. 


is the rotational diffusion coefficient for a single sphere rotating 
with Monte Carlo dynamics. In Eq. 17 0nis(Ar) is the mean- 


squared rotation angle (about a fixed arbitrary axis centered at 
the center of mass) that a sphere would experience if we applied 
our Stokes-scaled VMMC algorithm to volume elements within 
the sphere. This factor is necessary for rotations because, un¬ 
like for translations, the hydrodynamic radius (Eq. [n) depends 
on the center and axis of rotation. Defining 0(Ar,n,r) as the 
center-of-mass rotation angle for a rotation by Ar around n cen¬ 
tered at r, we find 


On 


s(\) - 2 ^ 


(6(Ar, ft, ^center)) 


Ro 


Ru{rc,n) 


n,\rc\<Ro 

(18) 


where 


Rnircn) = 


(—[ I 

V 27tRl Jr<RQ 


l(r-rc) X n\ 


1/2 


(19) 


is the center- and axis-dependent hydrodynamic radius. In the 
limit Ar 0, we numerically calculated ems(Ar)^e0s-0.14. 


As shown in Eig.[^(c), we find that for a range of tetrahedral 
cluster sizes Dr saturates from below near a value 
for Ar > 1. (Dr drops again for A^ > 5 due to the inability to 
resolve rotations when individual rotations exceed n radians.) 
We chose to parameterize Ar = 1 to ensure that rotations are 
saturated at the large-Ar limit, allowing translations to accom¬ 
modate internal relaxation. 

Third, we adjusted the attempt probability for translation and 
rotation, pi and pr = I — pu to correct for the numerical dis¬ 
crepancies between the predicted and measured diffusion coef¬ 
ficients for tetrahedral clusters (the factors 0.6 for translations 
and 1.5 for rotations discussed above). We achieved this by 
performing VMMC simulations of a tetrahedra of four spheres 
with various translation attempt probabilities pt and comparing 
the relative diffusion coefficients to the expected relationship 
(see Eq. 

A 4^ 2 




( 20 ) 


If the prefactors for the diffusion coefficients followed Eq. 15 
and Eq.[T^ Eq.[20| would be satisfied for 




Pt 


At 


= c 


fRoAr 


( 21 ) 


we 


Inserting the Ap ^ 0 limit, 0nis(Ar) — 0.14, into Eq. 21 
find c 0.31 and pt — 0.6. Instead, Eig. (d) shows that 
Pt = 0.7 (corresponding to c = 0.4) results in better agreement 
with Eq.[^ We therefore fix pi via Eq. [21] with c = 0.4. 

Einally, we showed that fixing the parameters as above re¬ 
sults in diffusion coefficients agreeing with the Stokes solutions 
(Eq.§ for tetrahedra of various sizes and with various interac¬ 
tion ranges. We fixed the translational step size At = 4A/?o to 
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be proportional to the interaction range to ensure that the ratio 
of single-particle to whole-cluster moves be consistent across 
interaction ranges. We fixed Ar = 1 and pi = 0.4(/?oA^/A^)^ as 
described above. Together, this parameterization fixes the time 
per Monte Carlo cycle, 

6 9 

^cycle ~ ^0? (22) 

where to is the natural Brownian time scale (Eq.[^. Fig.p^(e) 
and (f) show that Di and follow the Stokes solutions (Eq.|^ 
over a broad range of cluster sizes and interaction ranges, in 
stark contrast to single-particle Monte Carlo simulations, which 
follow the much more strongly size-dependent free draining so¬ 
lutions Dt ^ and ^ for small step sizes. 
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